clearvars -except HCR
path1='C:\Users\Tom�s\Documents\1. Regehr Lab\Data\HCR\OTRcrex Ai14 Slc6a5 Aldh1a3 12.7.20\Processing\4. position measurements\';
dir_to_search=[path1];
% pathtemp='C:\Users\Tom�s\Documents\1. Regehr Lab\Data\HCR\Registration and analysis\3.Channel colocalization\s1b2\';
path2='C:\Users\Tom�s\Documents\1. Regehr Lab\Data\HCR\OTRcrex Ai14 Slc6a5 Aldh1a3 12.7.20\Processing\1. Processed\';
path3='C:\Users\Tom�s\Documents\1. Regehr Lab\Data\HCR\'

% get the folder contents
d = dir(path1);
% remove all files (isdir property is 0)
dfolders = d([d(:).isdir]==1);
% remove '.' and '..'
dfolders = dfolders(~ismember({dfolders(:).name},{'.','..'}));

clear Cells ch1values ch2values ch3values cumch1Mat cumch2Mat cumch3Mat

for j=1:numel(dfolders)

    Cells=HCR(j).CC(:,[1:2]);


ch1=imread([path2 dfolders(j).name '\' 'MAX_tdTomato_.tif'],1);% Fiji saved
ch2=imread([path2 dfolders(j).name '\' 'Slc6a5tdTomSubMatlabout.tiff'],1);% Matlab out
ch3=imread([path2 dfolders(j).name '\'  'Aldh1a3Matlabout.tiff'],1); %Matlab out 

% figure
% imagesc(ch1) 
% caxis([0 300])
% hold on
% scatter(CCs(:,1),CCs(:,2),10,'filled','r')

ch1Mat=nan(121,121,size(Cells,1));
ch2Mat=nan(121,121,size(Cells,1));
ch3Mat=nan(121,121,size(Cells,1));
for i=2:size(Cells,1)
    ch1Mat(:,:,i)=ch1([round(Cells(i,2))-60:round(Cells(i,2))+60],round([round(Cells(i,1))-60:round(Cells(i,1))+60]));
    ch1MatR(:,:,i)=imresize(ch1Mat(:,:,i),[61 61],'bicubic'); %resize to match figure 2 analysis
    
    ch2Mat(:,:,i)=ch2([round(Cells(i,2))-60:round(Cells(i,2))+60],round([round(Cells(i,1))-60:round(Cells(i,1))+60]));
    ch2MatR(:,:,i)=imresize(ch2Mat(:,:,i),[61 61],'bicubic');%resize to match figure 2 analysis

    ch3Mat(:,:,i)=ch3([round(Cells(i,2))-60:round(Cells(i,2))+60],round([round(Cells(i,1))-60:round(Cells(i,1))+60]));
    ch3MatR(:,:,i)=imresize(ch3Mat(:,:,i),[61 61],'bicubic');%resize to match figure 2 analysis

end

if j==1
    cumch1Mat=ch1MatR;
    cumch2Mat=ch2MatR;
    cumch3Mat=ch3MatR;
else
end

% concatenate across slices
if ~(j==1)
    cumch1Mat=cat(3,cumch1Mat,ch1MatR);
    cumch2Mat=cat(3,cumch2Mat,ch2MatR);
    cumch3Mat=cat(3,cumch3Mat,ch3MatR);
else
end





end


% %% Display  mosaic of top cells (3 rows, 3 colums, total 9 cells)
% for idx=1:9
% randidx(idx)=round(rand(1)*size(cumch1Mat,3));
% end
% cell1=cumch1Mat(:,:,randidx(1))
% cell2=cumch1Mat(:,:,randidx(2))
% cell3=cumch1Mat(:,:,randidx(3))
% cell4=cumch1Mat(:,:,randidx(4))
% cell5=cumch1Mat(:,:,randidx(5))
% cell6=cumch1Mat(:,:,randidx(6))
% cell7=cumch1Mat(:,:,randidx(7))
% cell8=cumch1Mat(:,:,randidx(8))
% cell9=cumch1Mat(:,:,randidx(9))
% 
% 
% row1=cat(1,cat(1,cell1,cell2,cell3));
% row2=cat(1,cat(1,cumch1Mat(:,:,randidx(4)),cumch1Mat(:,:,randidx(5)),cumch1Mat(:,:,randidx(6))));
% row3=cat(1,cat(1,cumch1Mat(:,:,randidx(7)),cumch1Mat(:,:,randidx(8)),cumch1Mat(:,:,randidx(9))));
% 
% mosaic1=cat(2,(cat(2,row1,row2)),row3);
% 
% row1=cat(1,cat(1,cumch2Mat(:,:,randidx(1)),cumch2Mat(:,:,randidx(2)),cumch2Mat(:,:,randidx(3))));
% row2=cat(1,cat(1,cumch2Mat(:,:,randidx(4)),cumch2Mat(:,:,randidx(5)),cumch2Mat(:,:,randidx(6))));
% row3=cat(1,cat(1,cumch2Mat(:,:,randidx(7)),cumch2Mat(:,:,randidx(8)),cumch2Mat(:,:,randidx(9))));
% 
% mosaic2=cat(2,(cat(2,row1,row2)),row3);
% 
% row1=cat(1,cat(1,cumch3Mat(:,:,randidx(1)),cumch3Mat(:,:,randidx(2)),cumch3Mat(:,:,randidx(3))));
% row2=cat(1,cat(1,cumch3Mat(:,:,randidx(4)),cumch3Mat(:,:,randidx(5)),cumch3Mat(:,:,randidx(6))));
% row3=cat(1,cat(1,cumch3Mat(:,:,randidx(7)),cumch3Mat(:,:,randidx(8)),cumch3Mat(:,:,randidx(9))));
% 
% mosaic3=cat(2,(cat(2,row1,row2)),row3);
% % 
% figure;
% imagesc(mosaic1)
% colormap([cmap('red',100,0,50)]);
% caxis([200 20000])
% 
% figure
% imagesc(mosaic2)
% colormap([cmap('cyan',100,0,50)]);
% caxis([200 1500])
% 
% figure
% imagesc(mosaic3)
% colormap([cmap('magenta',100,0,50)]);
% caxis([200 1500])
% 
% % 
%  %%
% imwrite(uint16(mosaic1),[path3 'mosaic.tif'])
% imwrite(uint16(mosaic2),[path3 'mosaic.tif'],'WriteMode','append')
% imwrite(uint16(mosaic3),[path3 'mosaic.tif'],'WriteMode','append')

%%
% make nice figure 
%display all the cells stacked together to look at data distribution

radius=32;
ylimitprofiles=270;
figure
c1Ax=subaxis(2,3,1)
  imagesc(nanmean(cumch1Mat,3))
  caxis([0 15000])
    y1=colormap(c1Ax,[cmap('gray',100,0,0)]);

  axis square
  axis off
  subaxis(2,3,4)
%   plot(mean(mean(cumch1Mat,3),1))
meanGreen=nanmean(cumch1Mat,3);
greenHor=mean(meanGreen([29:33],:),1);
greenVert=mean(meanGreen(:,[29:33]),2);
greenProf=mean([greenVert';greenHor],1)
  hold on
  plot(greenProf,'k')
%   plot(mean(mean(cumch1Mat,3),2)) %average profiles in x and y to decide on boundary of circular mask
  
  ylim([0 ylimitprofiles*40/2])
%   vline([30-radius 30+radius],'--k')
  %draw axis for visual purpose
  vline(0,'k')
  hline(0,'k')
  axis square
  axis off
  
 c2Ax= subaxis(2,3,2)
  imagesc(nanmean(cumch2Mat,3))
  caxis([0 ylimitprofiles])
  y2=colormap(c2Ax,[cmap('cyan',100,0,50)]);
  axis square
  axis off
  subaxis(2,3,5)
%   plot(mean(mean(cumch2Mat,3),1))
%   hold on
%   plot(mean(mean(cumch2Mat,3),2)) %average profiles in x and y to decide on boundary of circular mask
meanCyan=nanmean(cumch2Mat,3);
CyanHor=nanmean(meanCyan([29:33],:),1);
CyanVert=nanmean(meanCyan(:,[29:33]),2);
CyanProf=nanmean([CyanVert';CyanHor],1)
  hold on
  plot(CyanProf,'k')
   ylim([0 ylimitprofiles/2])
%    vline([30-radius 30+radius],'--k')
   %draw axis for visual purpose
  vline(0,'k')
  hline(0,'k')
  axis square
  axis off
  
   c3Ax= subaxis(2,3,3)
  imagesc(nanmean(cumch3Mat,3))
  
  caxis([0 ylimitprofiles])
    y3=colormap(c3Ax,[cmap('magenta',100,0,50)]);

  axis square
  axis off
  subaxis(2,3,6)
%   plot(mean(mean(cumch3Mat,3),1))
%   hold on
%   plot(mean(mean(cumch3Mat,3),2)) %average profiles in x and y to decide on boundary of circular mask
meanMagenta=nanmean(cumch3Mat,3);
MagentaHor=mean(meanMagenta([29:33],:),1);
MagentaVert=mean(meanMagenta(:,[29:33]),2);
MagentaProf=mean([MagentaVert';MagentaHor],1)
  hold on
  plot(MagentaProf,'k')
   ylim([0 ylimitprofiles/2])
%    vline([30-radius 30+radius],'--k')
   %draw axis for visual purpose
  vline(0,'k')
  hline(0,'k')
  axis square
  axis off
 
  %% masking
    
% Create a logical image of a circle with specified
% diameter, center, and image size.
% First create the image.
imageSizeX = 61;
imageSizeY = 61;
[columnsInImage, rowsInImage] = meshgrid(1:imageSizeX, 1:imageSizeY);
% Next create the circle in the image.
centerX = ceil(imageSizeX / 2); % in the middle.
centerY = ceil(imageSizeY / 2);
radius = 50; %based on profiles

circlePixelsLogical = (rowsInImage - centerY).^2 ...
    + (columnsInImage - centerX).^2 <= radius.^2;
% circlePixels is a 2D "logical" array.
% Now, display it.
figure
image(circlePixelsLogical) ;
colormap([0 0 0; 1 1 1]);
title('Binary image of a circle');
axis square;

circlePixels01=double(circlePixelsLogical);
test=circlePixels01.*cumch1Mat; %mask out crap outisde cells 
figure
imagesc(nanmean(test,3))
caxis([0 50])
axis square



%% take averages of each cell for each channel inside the circular mask
clear ch1values
clear ch2values
clear ch3values
clear ch1tsne
clear ch2tsne
clear ch3tsne
for c=1:size(cumch1Mat,3)
    tempdata1=cumch1Mat(:,:,c);
%     ch1tsne(c,:)=tempdata1(circlePixelsLogical);
    ch1values(c)=mean(tempdata1(circlePixelsLogical));
    
     tempdata2=cumch2Mat(:,:,c);
%       ch2tsne(c,:)=tempdata2(circlePixelsLogical);
    ch2values(c)=mean(tempdata2(circlePixelsLogical));
    
     tempdata3=cumch3Mat(:,:,c);
%       ch3tsne(c,:)=tempdata3(circlePixelsLogical);
    ch3values(c)=mean(tempdata3(circlePixelsLogical));
    
end

figure
limVal=500;
subaxis(1,3,1)
scatter(ch1values,ch2values,10,'filled','k')
xlabel('tdTomato')
ylabel('Slc6a5')
axis square
xlim([0 30e3])
ylim([0 limVal])
subaxis(1,3,2)
scatter(ch1values,ch3values,10,'filled','k')
xlabel('tdTomato')
ylabel('Aldh1a3')
xlim([0 30e3])
ylim([0 limVal])
axis square
subaxis(1,3,3)
scatter(ch2values,ch3values,10,'filled','k')
xlabel('Slc6a5')
ylabel('Aldh1a3')
xlim([0 limVal])
ylim([0 limVal])
axis square

%% new section to try out the combination of circular mask with threshold mask based on tdTomato fluorescence


for i=1:size(cumch1Mat,3)
a=cumch1Mat(:,:,i);
back=prctile(a,5);
maxCelle=max(max(a([22:37],[22:37])));
b=a>back*3 & a<maxCelle;
threshch1Masks(:,:,i)=b & circlePixelsLogical;

% imshowpair(double(b),a,'montage')
% pause
end

% take mean values
for c=1:size(cumch1Mat,3)
    tempdata1=cumch1Mat(:,:,c);
%     ch1tsne(c,:)=tempdata1(circlePixelsLogical);
    ch1valuesT(c)=mean(tempdata1(threshch1Masks(:,:,c)));
    
     tempdata2=cumch2Mat(:,:,c);
%       ch2tsne(c,:)=tempdata2(circlePixelsLogical);
    ch2valuesT(c)=mean(tempdata2(threshch1Masks(:,:,c)));
    
     tempdata3=cumch3Mat(:,:,c);
%       ch3tsne(c,:)=tempdata3(circlePixelsLogical);
    ch3valuesT(c)=mean(tempdata3(threshch1Masks(:,:,c)));
    masksize(c)=sum(sum(threshch1Masks(:,:,c)));
    
end

figure
limVal=500;
tdTomlog=ch1values>2e3;
subaxis(1,3,1)
scatter(ch1valuesT(tdTomlog),ch2valuesT(tdTomlog),10,'filled','k','MarkerFaceAlpha',.7,'MarkerEdgeAlpha',.7)
xlabel('tdTomato')
ylabel('Slc6a5')
axis square
xlim([0 30e3])
ylim([0 limVal])
subaxis(1,3,2)
scatter(ch1valuesT(tdTomlog),ch3valuesT(tdTomlog),10,'filled','k','MarkerFaceAlpha',.7,'MarkerEdgeAlpha',.7)
xlabel('tdTomato')
ylabel('Aldh1a3')
xlim([0 30e3])
ylim([0 limVal])
axis square
subaxis(1,3,3)
scatter(ch2valuesT(tdTomlog),ch3valuesT(tdTomlog),10,'filled','k','MarkerFaceAlpha',.7,'MarkerEdgeAlpha',.7)
xlabel('Slc6a5')
ylabel('Aldh1a3')
xlim([0 limVal])
ylim([0 limVal])
axis square


%% Remake average HCR plot with cells above a certain tdTomato threshold

%%
% make nice figure 
%display all the cells stacked together to look at data distribution

tdTomlog=ch1values>2e3;
ylimitprofiles=300;
figure
c1Ax=subaxis(2,3,1)
  imagesc(nanmean(cumch1Mat(:,:,tdTomlog),3))
  caxis([0 20000])
    y1=colormap(c1Ax,[cmap('gray',100,0,0)]);

  axis square
  axis off
  subaxis(2,3,4)
%   plot(mean(mean(cumch1Mat,3),1))
meanGreen=nanmean(cumch1Mat(:,:,tdTomlog),3);
greenHor=nanmean(meanGreen([29:33],:),1);
greenVert=nanmean(meanGreen(:,[29:33]),2);
greenProf=nanmean([greenVert';greenHor],1)
  hold on
  plot(greenProf,'k')
%   plot(mean(mean(cumch1Mat,3),2)) %average profiles in x and y to decide on boundary of circular mask
  
%   ylim([0 ylimitprofiles*40/2])
%   vline([30-radius 30+radius],'--k')
  %draw axis for visual purpose
  ylim([0 20000])
  vline(0,'k')
  hline(0,'k')
  axis square
  axis off
  
 c2Ax= subaxis(2,3,2)
  imagesc(nanmean(cumch2Mat(:,:,tdTomlog),3))
  caxis([0 ylimitprofiles])
  y2=colormap(c2Ax,[cmap('cyan',100,0,50)]);
  axis square
  axis off
  subaxis(2,3,5)
%   plot(mean(mean(cumch2Mat,3),1))
%   hold on
%   plot(mean(mean(cumch2Mat,3),2)) %average profiles in x and y to decide on boundary of circular mask
meanCyan=nanmean(cumch2Mat(:,:,tdTomlog),3);
CyanHor=nanmean(meanCyan([29:33],:),1);
CyanVert=nanmean(meanCyan(:,[29:33]),2);
CyanProf=nanmean([CyanVert';CyanHor],1)
  hold on
  plot(CyanProf,'k')
   ylim([0 ylimitprofiles])
%    vline([30-radius 30+radius],'--k')
   %draw axis for visual purpose
  vline(0,'k')
  hline(0,'k')
  axis square
  axis off
  
   c3Ax= subaxis(2,3,3)
  imagesc(nanmean(cumch3Mat(:,:,tdTomlog),3))
  
  caxis([0 ylimitprofiles])
    y3=colormap(c3Ax,[cmap('magenta',100,0,50)]);

  axis square
  axis off
  subaxis(2,3,6)
%   plot(mean(mean(cumch3Mat,3),1))
%   hold on
%   plot(mean(mean(cumch3Mat,3),2)) %average profiles in x and y to decide on boundary of circular mask
meanMagenta=nanmean(cumch3Mat(:,:,tdTomlog),3);
MagentaHor=mean(meanMagenta([29:33],:),1);
MagentaVert=mean(meanMagenta(:,[29:33]),2);
MagentaProf=mean([MagentaVert';MagentaHor],1)
  hold on
  plot(MagentaProf,'k')
   ylim([0 ylimitprofiles])
%    vline([30-radius 30+radius],'--k')
   %draw axis for visual purpose
  vline(0,'k')
  hline(0,'k')
  axis square
  axis off




%% manually add values to cell types for cumulative plot (this is lazy but only needs tob e done once).

% 
% CCValues(:,1)=ch1values;
% CCValues(:,2)=ch2values;
% CCValues(:,3)=ch3values;
% 
% 
% GlobularValues(:,1)=ch1values;
% GlobularValues(:,2)=ch2values;
% GlobularValues(:,3)=ch3values;

% 
% 
% LugaroValues(:,1)=ch1values;
% LugaroValues(:,2)=ch2values;
% LugaroValues(:,3)=ch3values;


% 
% GolgiValues(:,1)=ch1values;
% GolgiValues(:,2)=ch2values;
% GolgiValues(:,3)=ch3values;


% MLI2Values(:,1)=ch1values;
% MLI2Values(:,2)=ch2values;
% MLI2Values(:,3)=ch3values;
% 


% %% Make cumulative plots
% MLI2Color=[233 83 30]./255; % MLI2 color
% CCColor=[30, 180, 233]./255; % CCs color from illustrator
% GlobularColor=[74 184 89]./255; %Globular color
% 
% LugaroColor=[127 63 152]./255; %from Illustrator
%  GolgiColor=[177 179 182]./255; %Golgi color illustrator
% 
% 
% 
% binWidth=2;
% figure
% hold on
% Ch1CChist=histogram(CCValues(:,1),'binwidth',binWidth,'FaceAlpha',0,'LineStyle', 'none','normalization' ,'cdf')
% Ch1Globularhist=histogram(GlobularValues(:,1),'binwidth',binWidth,'FaceAlpha',0,'LineStyle', 'none','normalization' ,'cdf')
% Ch1Lugarohist=histogram(LugaroValues(:,1),'binwidth',binWidth,'FaceAlpha',0,'LineStyle', 'none','normalization' ,'cdf')
% Ch1Golgihist=histogram(GolgiValues(:,1),'binwidth',binWidth,'FaceAlpha',0,'LineStyle', 'none','normalization' ,'cdf')
% Ch1MLI2hist=histogram(MLI2Values(:,1),'binwidth',binWidth,'FaceAlpha',0,'LineStyle', 'none','normalization' ,'cdf')
% % Globularhist=histogram(GlobularPCLdistance./GlobularMLlength,'binwidth',binWidth,'FaceAlpha',0.5,'LineStyle', 'none')
% 
% 
% 
% figure
% hold on
% Ch1CCStair=stairs(Ch1CChist.BinEdges(1:end-1), Ch1CChist.Values,'color',CCColor,'linewidth',2)
% Ch1GlobularStair=stairs(Ch1Globularhist.BinEdges(1:end-1), Ch1Globularhist.Values,'color',GlobularColor,'linewidth',2)
% Ch1LugaroStair=stairs(Ch1Lugarohist.BinEdges(1:end-1), Ch1Lugarohist.Values,'color',LugaroColor,'linewidth',2)
% Ch1GolgiStair=stairs(Ch1Golgihist.BinEdges(1:end-1), Ch1Golgihist.Values,'color',GolgiColor,'linewidth',2)
% Ch1MLI2Stair=stairs(Ch1MLI2hist.BinEdges(1:end-1), Ch1MLI2hist.Values,'color',MLI2Color,'linewidth',2)
% 
% 
% figure
% hold on
% Ch2CChist=histogram(CCValues(:,2),'binwidth',binWidth,'FaceAlpha',0,'LineStyle', 'none','normalization' ,'cdf')
% Ch2Globularhist=histogram(GlobularValues(:,2),'binwidth',binWidth,'FaceAlpha',0,'LineStyle', 'none','normalization' ,'cdf')
% Ch2Lugarohist=histogram(LugaroValues(:,2),'binwidth',binWidth,'FaceAlpha',0,'LineStyle', 'none','normalization' ,'cdf')
% Ch2Golgihist=histogram(GolgiValues(:,2),'binwidth',binWidth,'FaceAlpha',0,'LineStyle', 'none','normalization' ,'cdf')
% Ch2MLI2hist=histogram(MLI2Values(:,2),'binwidth',binWidth,'FaceAlpha',0,'LineStyle', 'none','normalization' ,'cdf')
% % Globularhist=histogram(GlobularPCLdistance./GlobularMLlength,'binwidth',binWidth,'FaceAlpha',0.5,'LineStyle', 'none')
% 
% 
% 
% figure
% hold on
% Ch2CCStair=stairs(Ch2CChist.BinEdges(1:end-1), Ch2CChist.Values,'color',CCColor,'linewidth',2)
% Ch2GlobularStair=stairs(Ch2Globularhist.BinEdges(1:end-1), Ch2Globularhist.Values,'color',GlobularColor,'linewidth',2)
% Ch2LugaroStair=stairs(Ch2Lugarohist.BinEdges(1:end-1), Ch2Lugarohist.Values,'color',LugaroColor,'linewidth',2)
% Ch2GolgiStair=stairs(Ch2Golgihist.BinEdges(1:end-1), Ch2Golgihist.Values,'color',GolgiColor,'linewidth',2)
% Ch2MLI2Stair=stairs(Ch2MLI2hist.BinEdges(1:end-1), Ch2MLI2hist.Values,'color',MLI2Color,'linewidth',2)
% 
% 
% figure
% hold on
% Ch3CChist=histogram(CCValues(:,3),'binwidth',binWidth,'FaceAlpha',0,'LineStyle', 'none','normalization' ,'cdf')
% Ch3Globularhist=histogram(GlobularValues(:,3),'binwidth',binWidth,'FaceAlpha',0,'LineStyle', 'none','normalization' ,'cdf')
% Ch3Lugarohist=histogram(LugaroValues(:,3),'binwidth',binWidth,'FaceAlpha',0,'LineStyle', 'none','normalization' ,'cdf')
% Ch3Golgihist=histogram(GolgiValues(:,3),'binwidth',binWidth,'FaceAlpha',0,'LineStyle', 'none','normalization' ,'cdf')
% Ch3MLI2hist=histogram(MLI2Values(:,3),'binwidth',binWidth,'FaceAlpha',0,'LineStyle', 'none','normalization' ,'cdf')
% % Globularhist=histogram(GlobularPCLdistance./GlobularMLlength,'binwidth',binWidth,'FaceAlpha',0.5,'LineStyle', 'none')
% 
% 
% 
% figure
% hold on
% Ch3CCStair=stairs(Ch3CChist.BinEdges(1:end-1), Ch3CChist.Values,'color',CCColor,'linewidth',2)
% Ch3GlobularStair=stairs(Ch3Globularhist.BinEdges(1:end-1), Ch3Globularhist.Values,'color',GlobularColor,'linewidth',2)
% Ch3LugaroStair=stairs(Ch3Lugarohist.BinEdges(1:end-1), Ch3Lugarohist.Values,'color',LugaroColor,'linewidth',2)
% Ch3GolgiStair=stairs(Ch3Golgihist.BinEdges(1:end-1), Ch3Golgihist.Values,'color',GolgiColor,'linewidth',2)
% Ch3MLI2Stair=stairs(Ch3MLI2hist.BinEdges(1:end-1), Ch3MLI2hist.Values,'color',MLI2Color,'linewidth',2)
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 

% CCch1=ch1tsne; 156
% CCch2=ch2tsne; 55
% CCch3=ch3tsne; 607

% MLI2tsne=horzcat(ch1tsne./156,ch2tsne./55,ch3tsne./607);
% Golgitsne=horzcat(ch1tsne./156,ch2tsne./55,ch3tsne./607);
% CCtnse=horzcat(CCch1./156,CCch2./55,CCch3./607);
% Globulartsne=horzcat(ch1tsne./156,ch2tsne./55,ch3tsne./607); %uses masked circle pixels
%% Use t-sne on the masked images with appropriate labels too see clustering of different cell types.
% alternatively use average value of masked image for each channel and
% normalized position in ML.

%make X vector for tsne, it will be size(X)= numberofCells by linearized
%pixelspercell
%group matrix g is going to be size = number of cells by 1
% X=vertcat(CCtsne,Globulartsne,Golgitsne,MLI2tsne);
% g=vertcat(ones(size(CCtsne,1),1),2*ones(size(Globulartsne,1),1),3*ones(size(Golgitsne,1),1),4*ones(size(MLI2tsne,1),1))
% rng default % for reproducibility
% Y = tsne(X,'Algorithm','barneshut','NumPCAComponents',0,'perplexity',500);
% figure;
% gscatter(Y(:,1),Y(:,2),g)
%  
    